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A variational calculation of the energy levels of the class of PT-invariant quantum mechani- 
cal models described by the non-Hermitian Hamiltonian H = p 2 — (ix) N with N positive and x 
complex is presented. The energy levels are determined by finding the stationary points of the 
functional (H)(a,b,c) = (f c dx tp{x)Hilj{x)^ / (J c dx ip 2 (x)j , where ip(x) — (ix) c exp (a(ix) b ) is a 
three-parameter class of "PT- invariant trial wave functions. The integration contour C used to define 
{H)(a,b, c) lies inside a wedge in the complex- £ plane in which the wave function falls off exponen- 
tially at infinity. Rather than having a local minimum the functional has a saddle point in the 
three-parameter (a, b, c)-space. At this saddle point the numerical prediction for the ground-state 
energy is extremely accurate for a wide range of N. The methods of supersymmetric quantum 
mechanics are used to determine approximate wave functions and energy eigenvalues of the excited 
states of this class of non-Hermitian Hamiltonians. 



I. INTRODUCTION 

In a recent letter jl| the spectra of the class of non-Hermitian "PT-symmetric Hamiltonians of the form 

H = p 2 - (ix) N {N>2) (1.1) 

were shown to be real and positive. It has been conjectured that the reality and positivity of the spectra are a 
consequence of VT symmetry. 

The Schrodinger differential equation corresponding to the eigenvalue problem Hijj = Etj) is 

- i)"{x) - (ix) N ip(x) = E^{x). (1.2) 

To obtain real eigenvalues from this equation it is necessary to define the boundary conditions properly. This was 
done in Ref. |l[ by analytically continuing in the parameter N away from the harmonic oscillator value N = 2. This 
analytic continuation defines the boundary conditions in the complex- a; plane. The regions in the cut complex- x 
plane in which tJj(x) vanishes exponentially as |x| — > oo are wedges. The wedges for N > 2 were chosen to be the 
continuations of the wedges for the harmonic oscillator, which are centered about the negative and positive real axes 
and have angular opening J. For arbitrary N > 2 the anti-Stokes' lines at the centers of the left and right wedges lie 
below the real axis at the angles 
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'right 
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N + 2 2 



(1.3) 



The opening angle of these wedges is ^ 2 . In Ref. Q the time- independent Schrodinger equation was integrated 
numerically inside the wedge to determine the eigenvalues to high precision. Observe that as N increases from its 
harmonic oscillator value (N = 2), the wedges bounding the integration path undergo a continuous deformation as a 
function of iV". 



In this work we show that the ground-state energy for the VT- invariant quantum mechanical models in Eq. (1.1) 
can be obtained from the condition that the functional 

, ,, x \n dx ib(x)Hib(x) 
E(a, b ,c,...) = { H)(aAc) = Jc Jl> 2 ^ ) , (1.4) 

be st atio nary as a function of the variational parameters a, b, c, ... of the trial ground-state wave function ip(x). In 
Eq. (1.4) the integration contour C lies in appropriate wedges in the complex plane, as explained in Sec. ||. 
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We obtain variational approximations to the higher eigenvalues and wave functions using the following procedure 
0,0: First, we obtain a 'PT-symmetric superpotential from the approximate ground-state trial wave function. From 
this we construct a supersymmetric 'PT-symmetric partner potential. Next, we use variational methods to find the 
ground-state energy and wave function of the Hamiltonian associated with the partner potential and, from that, a 
second superpotential. Iterating this process, we determine the higher-energy eigenvalues and the associated excited 
states of the original Hamiltonian. 

We now review this procedure in more detail. Let us consider the Hamiltonian whose eigenvalues E n °^ and 
eigenfunctions ipn^ are labeled by the index n. Subtracting the ground-state energy Eq ^ from the Hamiltonian 
gives a new Hamiltonian = — Eq°\ The eigenvalues of this new Hamiltonian are shifted accordingly, 



£q , but the eigenfunctions of remain unchanged: tpn (x) = ip n (x) 



E W - r(°) 

By construction, ifW has zero ground-state energy. Thus, we may regard as the first_component of a two- 
component supersymmetric PT-symmetric Hamiltonian that can be written in factored form 



H"> = A™ A™ 

-[IS-"^W1I7S + "^W1- 

The ground-state wave function of H W satisfies the first-order differential equation 

[~ + W ( - 1 \x)}¥ 1) (x)=0. 



From this equation we determine the superpotential W^ 1 



W (1 \x) = 



(0)'/ 



(x) a,™ (x) 



(0), 



(1.5) 



(1.6) 



(1.7) 



The supersymmetric partner of which we call H^> 

opposite order from that in Eq. (1.5): 



is given by the product of the operators A± in the 



_d?_ 

dx 2 



W^{x) 



W w (x) 



(1.8) 



Note that the spectrum of is the same as the spectrum of H^- 1 ' except that it lacks a state corresponding to 
the ground state of It is easy to verify that the eigenstates of H^> are explicitly given by ipi^ — A^'ip^}^ an d 



(2) iii 

the corresponding eigenvalues are En = E^! x . Thus, if we can find the ground-state wave function and energy of 
H( 2 \ then we can use these quantities to express the first excited wave function and energy of H (°) : 

(0) A^4 2 \x) 
n 0*0 = . 



E, 



(2) 



E (0) = E (2) 



E, 



(0) 



(1.9) 



By iterating this procedure we can construct a hierarchy of Hamiltonians H^" 1 ' (m =1,2,3, . . .). The successive 
energy levels of iJ (0) are then obtained by finding the ground-state wave functions and energies of if' ' using 
variational methods that we develop in Sec. |l]. It is important to realize that because the higher wave functions of 

H (0) 

require an increasing number of derivatives, we need increasingly accurate ground-state wave functions to obtain 
a reasonable approximation to the higher eigenvalues. A more extensive discussion of this procedure can be found in 
Refs.®. 

This paper is organized very simply. In Sec, fl we use variational methods to calculate the ground-state energy and 
wave function of the Hamiltonian H in Eq. (1.1) for various values of N. Then in Sec. EH we calculate the first few 
excited energy levels and wave functions. In Sec. ^ we make some concluding remarks and suggest some areas for 
future investigation. 
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II. VARIATIONAL ANSATZ FOR PT-SYMMETRIC HAMILTONIANS 



In previous variational calculations || two-parameter "post-Gaussian" nodeless trial wave functions of the general 
form 



ip(x) = exp (-a\x\ b ) (2.1) 

were used to obtain estimates for the ground states of the hierarchy of Hamiltonians constructed from the anhar- 
monic oscillator p 2 + x 4 . In this paper we obtain estimates for the ground states of the hierarchy of PT-symmetric 



Hamiltonians arising from H in Eq. (1.1). 

The concept of VT symmetry is discussed in detail in Refs. 0] and ||. To be precise, a function F(x) of a complex 
argument x is VT symmetric if 

[F(x)]*=F(-x*). (2.2) 
Thus, any real function of ix is manifestly VT symmetric. In previous numerical calculations |l| it was found that 



the eigenfunctions of the 'PT-symmetric Hamiltonians in Eq. (1.1) were all VT symmetric. Thus, for our variational 



calculation we choose a three-parameter trial wave function that is explicitly VT symmetric: 

i>{x) = {ix) c exp (a(ix) b ) , (2.3) 

where the variational parameters a, b, and c are real. 

The most notable difference between the Hermitian trial wave function in Eq. J2.l|) and the PT-symmetric trial 



wave fun ction in Eq. (2.3) is the number of variational parameters. It is not advantageous to have a prefactor in 
Eq. (2.1) of the form \x\ c because the expectation value of the kinetic term —d 2 /dx 2 of the Hamiltonian does not 



exist unless c > h. Because of the singularity at x — we must exclude the case c < \ in conventional variational 
calculations. However, positive values of c are in conflict with the large-|a;| behavior of the WKB approximation to 
the exact wave function; WKB theory predicts a negative value for the parameter c. 

We emphasize str ongl y that in conventional variational calculations the integration path is the real axis. For the 
case of i\){x) in Eq. (fyj) one cannot deform the contour of integration to avoid the origin because the function \x\ is 
not analytic. More generally, the notion of path independence is not applicable in conventional variational calculus 
because the functional to be minimized involves an integral over the absolute square tl>*(x)ifi(x) of the trial wave 



function. For the case of "PT-symmetric quantum mechanics the trial wave function ip(x) in Eq. (2.3) is analytic in 
the cut x plane. (The cut is taken to run along the positive-imaginary axis; this choice is required by VT symmetry.) 
Moreover, the functional involves an integral over the square, not the absolute square, of the wave function ||. These 
two properties imply that the integration contour can be deformed so long as the endpoints of the contour lie inside 



of the wedge in Eq. (1.3). In particul ar, the contour may be chosen to avoid the infinite singularity at the origin that 
occurs when the parameter c in Eq. (p^q) is negative. Thus, in PT-symmetric quantum mechanics we can include the 
additional parameter c anticipating that we will obtain highly accurate numerical results for the energy levels. 

In general, for Hermitian Hamiltonians it is well known that the variational method always gives an upper bound for 
the ground-state energy. Th at is , the stationary point of the expectation value of the Hamiltonian in the normalized 
trial wave function [see Eq. ( |l.4| )] is an absolute minimum. By contrast, for PT-invariant quantum mechanics, too 
little is known about the state space to prove any theorems about the nature of stationary points. Thus, we will 



simply look for all stationary points of the function E(a, b, c) in Eq. (tL4J), and not jus t minima 



The first step in our calculation is to evaluate the functional E(a,b,c) in Eq. (1.4). To do this we must make a 
definite choice of contour C . It is convenient to c hoose the path of integration in the complex- x plane to follow the 
rays for which the trial wave function in Eq. (2JJ) does not oscillate and falls off exponentially as exp (— ar b ). (We 



assume implicitly here that a is positive.) These rays are symmetrically placed about the negative- imaginary axis. 
Specifically, we set 

z = re idL , where 6 L = -- - - (2.4) 

2 b J 

on the left side of the imaginary axis, with the path running from r = -oo to r = (9l fixed), and 

z = re l0R , where 6 R = -| + | (2.5) 

on the right side of the imaginary axis, with the path running from r = to r — oo (Or fixed). Note that when 6 = 2, 
the path lies on the real axis; for b ^ 2, the path develops an elbow at the origin. 
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Our key assumption is that this integration path lie inside the asymptotic wedges described in Eq. (1.3). This 
requirement restricts the value of b: 

-(N + 2) < b < N + 2. (2.6) 
3 

Now, we evaluate the integral in Eq. fll>4Q in terms of Gamma functions: 

(p-f3 2 ~ 7 )r(i -0- iWP 



4/3 2 r(l + f3 - 7) 
r(l - f3 - j)a~ NI3 



(2.7) 



T{l-P-Np-i) 
where 

a = 2a, (3 = 1/6, 7 = 2c/b. (2.8) 

To find a stationary point of the function E(a, (3, 7) we must calculate its partial derivatives with respect to the 
parameters a, (3, and 7, and determine the values of these parameters for which the partial derivatives vanish. It is 
simplest to calculate the derivative with respect to a; requiring that this derivative be zero gives an expression for a 
in terms of (3 and 7: 

a (2 +N) p = 2^r(i + /?- 7 ) , 9) 

^-f3 + 1 )T{l-(3-N(3- 1 y [ - } 

Recalling that the parameter a, and therefore a, must be positive gives regions in the (/?, 7) plane of the allowable 
values of (3 and 7. These regions are shown as shaded areas in Fig. [l]. 

There are many distinct regions for which a is positive. We describe the boundaries of these regions in detail. First, 
all of these regions must satisfy the inequality 

' <P<^~o, (2-10) 



N + 2 -' ~N + 2' 

which follows from Eq. (2.6). These bounds are shown as two heavy vertical lines in Fig. [jj 
Second, there is an inverted parabola 

7-/5-/J 2 . (2.11) 

The lowest region is bounded above by this parabola and the next higher region is bounded below by this parabola. 
Third, there is an infinite sequence of downward-sloping straight lines of the general form 

j=k-/3(N + l) (k = 1, 2, 3, ...). (2.12) 

These lines bound the shaded regions on the left and on the right. 

Fourth, there is an infinite sequence of upward-sloping straight lines of the general form 

7 = fc + /3 (k = 1, 2, 3, ...). (2.13) 

These lines are upper and lower bounds of the shaded regions. 
WKB theory predicts that 

A 2 N 

(3= , 7 = . (2.14) 

1KT 1 O ' 1\.T 1 O v ' 



N + 2 1 ' N + 2 1 ' N + 2 

The special WKB values ((3 = jyq^,7 = —jrp^) are indicated by a heavy dot in Fig. [l[ Note that this point lies on 
the boundary of the only shaded region that lies below the (3 axis. We will restrict our attention to the stationary 
point (there is only one) of E(a,f3,j) that lies in this special region because it is nearest to the WKB point. We 
will see that the stationary point in this special region is a local maximum as a function of (3 and 7, and not a local 
minimum as in the case of conventional Hermitian theories. 

To locate the stationary point in the (j3, 7) plane we solve Eq. (2.9) for a and substitute the result back into the 
energy functional E(a.,(3,"f) in Eq. fl2.7|). We then construct a highly-accurate contour plot of the resulting function 
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of (3 and 7. From our analysis we are able to locate the stationary points to four significant digits for various values 
of TV. Our results for TV — 3, TV — 4, and TV = 5 are shown in Tables | and |n[ In Table III we show that in the (/?, 7) 
plane, the stationary point is a local maximum. 

To verify the accuracy of t he ground-state trial wave function we have calculated the normalized expectation value 
x p using ip(x) in Eq. ( |2.3| ): 

P»/,2/ 



Of -p 



x p |0) _ f c dxx p ip 2 (x) 



(0|0> 



J c ip 2 (x) 

p r(l-/3- 7 )q-^ 
lj r(l-/3-P/3- 7 )' 



(2.15) 



3, TV = 4, and TV = 5 we display in Tabic [V the 



Substituting the values of a, 0, and 7 taken from Table Q for TV 

expectation values of x p for P = 1, 2, . . . , 5. To measure the accuracy of our variational wave function we can check 

in two ways. First, we have calculated the exact values of (0| x |0)/ (0|0) for TV = 3, TV = 4, 
The values are — 0.590073i, — 0.866858i, and — 1.013102z, respectively. These numbers differ from those 
in the first row in Tabic ^ by about one part in a thousand. Second, combining the operator equation of motion 
obtained from H in Eq. ( Ll| ) , 



the numbers in Table IV 
and TV = 5 



x"(t) = 2iN(ix) 



N-l 



(2.16) 



with the time-translation invariance of (0| x |0), we find that the expectation value of x N 1 vanishes for an (ix) N 
theory. This is evident to a high degree of accuracy in Table |v|. 

One further check of the accuracy of the trial wave f uncti on is to compare its asymptotic behavior with that of the 
WKB approximation to il>(x). For example, from Eq. ( 2. 14f ) the WKB value of /3 for TV = 3 is 0.4. For a simpler two- 
parameter trial wave function in which we set 7 = 0, the stationary point occurs at (3 = 0.36. For the three-parameter 
case we get j3 = 0.3855, whic h is much closer to the value 0.4. Indeed, all of the values in Table § are very close to 



the WKB values in Eq. (2.14) 



III. HIGHER ENERGY LEVELS 

We follow the procedure outlined in Sec. | and use trial wave functions of the form in Eq. (|1|) to construct the first 
few excited states. As one might have expected, the numerical accuracy decreases for the higher states. The ground- 
state values of (a, j3, 7) for the case TV = 3 are (0.72066, 0.3855, —0.157), as shown in Table |. For the first excited state 
we obtain the values (0.55750,0.3568,-0.329) and for the second excited state we obtain (0.39925,0.3255,-0.466). 
Using these parameters we predict that the first excited energy level E\ is 4.11738 and that the second excited energy 
level E2 is 7.53886. The exact values of E\ and E2 obtained by numerical integration of the Schrodinger equation 



(1.2) (see Rcf. |1[) arc 4.1092 and 7.5621. Thus, the relative errors in our calculation of Eq, E%, and E% are roughly 
1 part in 2000, 4 parts in 2000, and 7 parts in 2000. 

For the case TV = 4 the ground-state values of {a,f3, 7) are (0.5800,0.3185,-0.2745) (see Table fl). For the first 
excited state we obtain (0.41763, 0.2915, —0.614). From these parameters we predict that the first excited energy level 
Ei is 6.03769. The exact value of Ei is 6.0033. Thus, the relative error in this calculation of E and Ei is roughly 3 
parts in 5000 and 25 parts in 5000. 



TABLE I. Variational calculation of the ground-state energy for the 'PT-symmetric Hamilto nian in Eq. (1.1) for three values 



of TV. The stationary point of E(a, /3, 7) for the three-parameter trial wave function in Eq. (2.3) is given for TV = 3, TV 
and TV = 5. 



N 


a 




7 


3 


0.7207 


0.3855 


-0.1570 


4 


0.5800 


0.3185 


-0.2745 


5 


0.4895 


0.2720 


-0.3610 



5 



TABLE II. Variational calculation of the ground-state energy for the 'PT-symmetric Hamiltonian in Eq. ( p_,l| ) for three values 
of N. Using the stationary point of E(a, (3, 7) taken from Table jjj we calculate the variational value of the ground-state e nerg y 
E va _ r and compare it with the exact value -Ecxact obtained by direct numerical integration of the Schrodinger equation in (1.2). 
Note that the relative error for N = 3, N = 4, and N = 5 is approximately one part in 2000. 

N E mT -Boxact rel. error 

3 1.156754 1.156267 0.042% 

4 1.478023 1.477149 0.062% 

5 1.909382 1.908265 0.059% 



TABLE III. The stationary point of E[a(/3, 7), j3, 7] is a local maximum in the (/?, 7) plane. To verify this the value of E for 

N = 4 is given at the stationary point j3 = 0.3185, 7 = —0.2745 and at four nearby points. In a conventional Hermitian theory 

the stationary point would be a local minimum. 

\ 7 E{a,p,~ 

0.3185 -0.2745 1.478023 

0.3185 -0.2750 1.478022 

0.3185 -0.2740 1.478022 

0.3186 -0.2745 1.478022 

0.3184 -0.2745 1.478022 



TABLE IV. Normalized expectation values of x p in the ground state. These expectation values were calculated from the 
trial wave function tp(x) in Eq. (Bjh with a, (3, and 7 given in Table fl] for iV = 3, iV = 4, and N = 5. 



p 


<0| x p \0) 




<0| z P |0) 




<0| x p |0> 






(010) 


JV=3 


(010) 


JV=4 


(010) 


N=5 



1 -0.590686i -0.867264i -1.013266i 

2 -0.000767 -0.518220 -0.864749 

3 -0.462705i 0.000869i 0.518217i 

4 -0.385363 -0.492631 0.002080 

5 -0.363227i 0.635997i 0.545538i 



Having verified that the variational calculation of the energies of the excited states is extremely accurate, we turn 
to the wave functions. We have performed two tests of the accuracy of the variational wave functions. First, we have 
computed the integral of the square of the wave functions. In general, we expect that the integral of the square of 
the nth excited wave function to be a real number having the sign pattern (—1)™. This is precisely what we find for 
the N — 3 wave functions (ipo, ipi, and ^2) and the N = 4 wave functions (ipo and ipi) that we have studied. To 
understand the origin of this alternating sign pattern recall the form of the wave functions for the harmonic oscillator 
[the N = 2 case of Eq. (|l .]])]. The harmonic oscillator is VT symmetric and all of its eigenfunctions can be made VT 
symmetric upon multiplying by the appropriate factor of i. The first three eigenfunctions in explicit 'PT-symmetric 
form are, apart from a real multiplicative constant, 

lH>{x) = e^> 2 / 2 , 
Mx) = («)e (,l)2/2 , 
Mx) = [2{ixf + l]e {lxS > 2 / 2 . 

These wave functions exhibit precisely this (—1)™ behavior. 

Second, we have found the nodes of the variational wave functions. For N = 3 the ground-state variational wave 
function is nodeless and the first excited state variational wave function has one node at x = — 0.703i. These results 
are extremely accurate: The exact ground-state wave function is also nodeless and the exact first-excited-state wave 
function has one node at x — — 0.533i. We observe the same qualitative features for the case N — 4; here, the node 
of the first excited state, as determined by our variational calculation, is located at x = — 0.972i. 



IV. CONCLUSIONS 



The results of this investigation suggest several intriguing avenues of research. The most surprising aspect of this 
variational calculation is that even though the stationary point is a saddle point in parameter space, and not a 
minimum, this saddle point clearly gives correct physical information. Clearly, it is the nonpositivity of the norm that 
accounts for the fact that the stationary point is a saddle point. The norm used in this paper is 

/ dxip 2 (x) (4.1) 
Jc 

rather than the conventional choice 

/ dxil)*{x)il>(x). (4.2) 
Jc 

However, the accuracy of the variational approach suggests that we have found an extremely powerful technique for 
obtaining approximate solutions to complex Sturm-Liouville eigenvalue problems. 

A possible explanation for the accuracy of our approach is that while the norm we are using is not positive definite 
in coordinate space, we believe that it is positive definite in momentum space. The positivity for our choice of norm is 
easy to verify in the special case of the harmonic oscillator. We believe that this positivity property does not undergo 
a sudden transition as we analytically continue the theory into the complex plane by increasing N. We point out 
that the strong advantage of our choice of norm is that ip (x) is an analytic function while ip* (x)ip(x) is not. This 
allows us to analytically continue Sturm-Liouville problems into the complex plane. (We point out that until now, 
Sturm-Liouville problems have been limited to the real axis.) 

An extremely interesting question to be considered in the future is the completeness of the eigenfunctions of a 
complex PT-symmetric Sturm-Liouville problem. It is easy to show that the eigenfunctions corresponding to different 
eigenvalues are orthogonal: 

dxil) m {x)^ n (x) — (m^n). (4.3) 

However, it is not known whether the eigenfunctions form a complete set and what norm should be used in this 
context. It is clear that in order to answer questions regarding completeness, a full understanding of the distribution 
of zeros of the wave functions must be achieved. For a Sturm-Liouville problem limited to the real axis completeness 
depends on the zeros becoming dense. For the case of a regular Sturm-Liouville problem the zeros interlace as they 
become dense. Our preliminary numerical investigations suggest that the zeros of the eigenfunctions of complex 
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PT-symmetric Sturm-Liouville problems also become dense in the complex plane and exhibit a two-dimensional 
generalization of the interlacing phenomenon. 

One area that we have investigated in great detail is the extension of these quantum mechanical variational calcula- 
tions to quantum field theory ||. The natural technique to use is based on the Schwinger-Dyson equations; successive 
truncation of these equations is in exact analogy with enlarging the space of parameters in the trial wave functions. 
We find that the Schwinger-Dyson equations provide an extremely accurate approximation scheme for calculating the 
Green's functions and masses of a PT-symmetric quantum field theory. 

Finally, we remark that the notion of probability in quantum mechanics is connected with the positivity of the 
norm in Sturm-Liouville problems. We note that the conventional quantum mechanical norm is neither an a naly tic 
function nor invariant with respect to PT symmetry, while the momentum-space version of the norm in Eq. (4.1) is 
both. It would be a remarkable advance if one could formulate a fully PT-symmetric quantum mechanics. 

We thank Dr. S. Boettcher for assistance in doing computer calculations. This work was supported in part by the 
U.S. Department of Energy. 
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FIG. 1. Regions (sh own as shaded) of acceptable values of /3 and 7. These values of /3 and 7 are acceptable because they 
give positive a in Eq. (E^J). 




